setwd("C:/Users/Kadeem Gilbert/Desktop/Microbiome project/ANCOM") # Set working directory
if (!require("vegan")) {install.packages("vegan"); require("vegan")}
# packages required for ANCOM
if (!require("ancom.R")) {install.packages("ancom.R"); require("ancom.R")}
if (!require("doParallel")) {install.packages("doParallel"); require("doParallel")}
if (!require("DT")) {install.packages("DT"); require("DT")}
if (!require("exactRankTests")) {install.packages("exactRankTests"); require("exactRankTests")}
if (!require("foreach")) {install.packages("foreach"); require("foreach")}
if (!require("ggplot2")) {install.packages("ggplot2"); require("ggplot2")}
if (!require("Rcpp")) {install.packages("Rcpp"); require("Rcpp")}
if (!require("shiny")) {install.packages("shiny"); require("shiny")}

# initializing otu data
botu.10 <- read.table("botu_h.txt",sep="\t",head=T,row.names=1) #read in OTU table
botu.notax <- botu.10[,2:ncol(botu.10)] # remove taxonomy
summary(rowSums(botu.notax))
summary(colSums(botu.notax))
botu.t <- t(botu.notax) # transpose rows and columns to get OTUs as columns, samples as rows
summary(rowSums(botu.t))
summary(colSums(botu.t))
botu.t <- botu.t[order(row.names(botu.t)),] #order row names
botu.t1k <- botu.t[rowSums(botu.t) > 1000,] # keep only samples with at least 1k seqs
botu.t1k.no0 <- botu.t1k[,colSums(botu.t1k) > 100] # keep only OTUs with nonzero counts above 100. 

# initializing metadata
meta <- read.table("MicrobiomeMetadata_ANCOM.txt",sep="\t",head=T,row.names=1) #read in OTU table
meta <- meta[order(row.names(meta)),] # order row names
meta.botu.t1k.no0 <- subset(meta, row.names(meta) %in% row.names(botu.t1k.no0)) #subsetting metadata to match filtered botus
row.names(meta.botu.t1k.no0) == row.names(botu.t1k.no0) # checking

# formatting data in "wide format"
Group <- meta.botu.t1k.no0$Morph
ID <- meta.botu.t1k.no0$SampleID
botu.t1k.gr.no0 <- cbind(botu.t1k.no0, Group)
botu.t1k.gr.no0 <- cbind(botu.t1k.gr.no0, ID)
botu.t1k.gr.no0 <- data.frame(botu.t1k.gr.no0)


write.table(botu.t1k.gr.no0, file = "botu.t1k.gr.nz.txt", sep = "\t", col.names = colnames(botu.t1k.gr.no0), row.names=FALSE, quote=FALSE)

# read it back in
# my.df2 <- read.table(file = "potu.t1k.gr.nz.txt",sep = ",", header = TRUE, stringsAsFactors = FALSE) 

# Running with the shiny application
# From there I follow upload instructions on the README pdf pages 5-6
shiny_ancom()

# Running ANCOM manually (README pdf pages 7-8)
ancom.p.out <- ANCOM( potu.t1k.gr.no01, sig = 0.05, multcorr = 2, repeated=TRUE )
ancom.p.out$W
ancom.p.out$detected
plot_ancom(ancom.p.out)

ancom.potu.dat <- read.table("Selected_OTU_data_potu.t1k.gr.nz.txt",sep="\t",head=T,row.names=1) #read in OTU table

#####################################################################################

# initializing otu data
eotu.10 <- read.table("eotu_h.txt",sep="\t",head=T,row.names=1) #read in OTU table
eotu.notax <- eotu.10[,2:ncol(eotu.10)] # remove taxonomy
summary(rowSums(eotu.notax))
summary(colSums(eotu.notax))
eotu.t <- t(eotu.notax) # transpose rows and columns to get OTUs as columns, samples as rows
summary(rowSums(eotu.t))
summary(colSums(eotu.t))
eotu.t <- eotu.t[order(row.names(eotu.t)),] #order row names
eotu.t1k <- eotu.t[rowSums(eotu.t) > 1000,] # keep only samples with at least 1k seqs
eotu.t1k.no0 <- eotu.t1k[,colSums(eotu.t1k) > 100] # keep only OTUs with nonzero counts above 100. If I don't do this step, the file is too big for the package

# initializing metadata
meta <- read.table("MicrobiomeMetadata_ANCOM.txt",sep="\t",head=T,row.names=1) #read in OTU table
meta <- meta[order(row.names(meta)),] # order row names
meta.eotu.t1k.no0 <- subset(meta, row.names(meta) %in% row.names(eotu.t1k.no0)) #subsetting metadata to match filtered eotus
row.names(meta.eotu.t1k.no0) == row.names(eotu.t1k.no0) # checking

# formatting data in "wide format"
Group <- meta.eotu.t1k.no0$pH
ID <- meta.eotu.t1k.no0$SampleID
eotu.t1k.gr.no0 <- cbind(eotu.t1k.no0, Group)
eotu.t1k.gr.no0 <- cbind(eotu.t1k.gr.no0, ID)
#rownames(eotu.t1k.gr.no0) <- c()
eotu.t1k.gr.no0 <- data.frame(eotu.t1k.gr.no0)

write.table(eotu.t1k.gr.no0, file = "eotu.t1k.gr.nz.txt", sep = "\t", col.names = colnames(eotu.t1k.gr.no0), row.names=FALSE, quote=FALSE)

# read it back in
# my.df2 <- read.table(file = "eotu.t1k.gr.nz.txt",sep = ",", header = TRUE, stringsAsFactors = FALSE) 

# Running with the shiny application
# From there I follow upload instructions on the README pdf pages 5-6
shiny_ancom()

# Running ANCOM manually (README pdf pages 7-8)
ancom.p.out <- ANCOM( eotu.t1k.gr.no01, sig = 0.05, multcorr = 2, repeated=TRUE )
ancom.p.out$W
ancom.p.out$detected
plot_ancom(ancom.p.out)

ancom.eotu.dat <- read.table("Selected_OTU_data_eotu.t1k.gr.nz.txt",sep="\t",head=T,row.names=1) #read in OTU table
